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1 .   INTRODUCTION 

Predictor-corrector  methods   furnish  attractive   algorithms 
for  the  solution  of  ordinary  differential  equations  because  of  the 
relatively  small  number  of  derivative  evaluations   required  per  step. 
For  example,   fourth  degree  predictor-corrector  methods   require  two 
derivative  evaluations  per  integration   step  while  the   corresponding 
Runge-Kutta  fourth   degree   algorithm  requires    four  derivative   evalu- 
ations per  integration  step.      The   difficulty  of  missing  starting 
values   is   typical   for  all  multistep  methods.      One   step  methods   such 
as  Runge-Kutta  methods   or  Taylor's   expansion   are   frequently  used  to 
find  missing  starting  values.      Other  starting  methods  which   are  more 
within  the   spirit   of  difference  methods    are   also  known   (Collatz,    "The 
Numerical  Treatment   of  Differential  Equations,"    3rd  Ed.    Springer, 
Berlin    (i960),   p.    8l). 

One   of  the  key  factors   to  be   considered  in  the  selection  of 
a  particular  predictor-corrector  pair  is  the  stability  of  the  numerical 
algorithm.      This   is  particularly   crucial  when  the   differential  equations 
being  solved  correspond  to   a  system  with   a  forcing  function  whose  time 
duration  or  period  is   relatively  long  compared  to  the  transient  time 
constants   of  the   system. 

Instability   is    also  introduced  by  the   fact  that   a  first  order 
equation  is   solved  with  the  help  of  higher  order  difference  equations. 
The  principal  root   of  the   difference   equation   approximates   the   true 
solution  of  the  equation.      Any  other  root   is   extraneous.      Aside   from 
these  parasitic   roots,  presence  of  nonlinearity  in  a  differential 
equation   can  cause   instability,  while  establishing  the   region  of  stability 
for  a  particular  method  this   is   generally  ignored. 


In  this  work  we  will  try  to  study  the  effect  of  nonlinearity 
in  an  ordinary  differential  equation.     Adams -Bashforth  predictor  and 
Adams-Moulton   corrector  of  orders  three  through  six  will  he  studied 
and  regions   of  their  absolute   and  relative   stability  established  for 
different  values   of  parameters.      This   study  is  not  exhaustive  because 
only  a  few  values   of  the  parameters  will  be   considered.      We  will 
introduce  the   following  terminology  to  describe  the  purpose  of  this 
study  in  more   detail. 

y'    =   f(x,  y)  a   <  x  <  b 

y(a)  =  y0 

g(x)   =   fy(x,  y(x)) 

gn  =   g(xn) 
hgn  =  h 

and  g,(xn)   =  <*g^ 

We  will  ignore  higher  derivatives   of  g(x).      Define   a  function   0  such  that 

0(o,  h)   =    [0  if  stable 
j 
•   1  if  unstable 

We  will  establish 

i.      The  region  of  numerical  stability  for  Adams -Bashforth 
predictor  and  Adams -Moult on   corrector  of  orders  three 
through  six  for  a  ^   0   and  compare   it  with  the   corresponding 
region  of  stability   for  a  =  0. 
ii.      It  was   conjectured  that   region  of  stability   for  a  ^  0 
for  a  particular  predictor-corrector  pair  lies   inside 
the   region  of  stability  for  a  =   0.      This  was   not   so  in 
some  examples . 
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iii.  Another  conjecture  we  wished  to  test  was  whether 

0(0,  h(l-J2h))  =  0 
for  j  =  0,  1,...,  k  for  a  k  step  method. 

=s>   0(a,  h)  =  0 
This  was  also  false. 
In  the  following  pages  these  three  properties  will  be 
studied  in  some  detail. 


2.   DESCRIPTION  OP  METHODS  USED  AND  RESULTS 

All  the  methods  we  use  for  solving  differential  equations 
have  various  possible  sources  of  error.   Circumstances  make  it  likely 
that  the  true  solution  differs  to  some  extent  from  our  computed 
solution,  and  we  try  to  make  the  discrepancy  as  small  as  required  for 
any  particular  application.   When  we  try  to  estimate  the  errors,  a 
number  of  factors  must  be  taken  into  account.   They  are 
i.   Inherent  instability  in  the  problem. 
ii.   Instability  introduced  by  our 
particular  method  of  solution, 
iii.   Complete  or  partial  neglect  of  the  truncation  error, 
iv.   The  presence  of  some  form  of  singularity, 
v.   Nonlinearity  in  a  differential  equation. 
The  problem  of  inherent  instability  is  independent  of  the 
method  of  solution.   We  will  not  consider  this  problem  in  this  work. 
We  will  mainly  be  concerned  with  the  effect  of  nonlinearity  on  the 
stability  of  predictor-corrector  methods.   To  get  trustworthy  results 
stability  with  a  large  margin  of  safety  under  all  circumstances  is 
essential.   In  lengthy  automatic  computations  the  number  of  elementary 
steps  may  be  as  large  as  10  or  more,  and  disturbances  in  unstable 
methods  typically  grow  exponentially  with  the  number  of  steps. 
Predictor-corrector  methods  may  cause  more  catastrophic  strong 
unstability  through  the  introduction  of  spurious  solutions .   A  result 
of  this  kind  is  hardly  surprising  since  the  higher  order  difference 
equation  must  introduce  spurious  solutions  having  no  relation  to  that 
of  the  first  order  differential  equation.   The  important  feature  is 
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that  the   solution  is  not  improved  by  a  reduction  in  the  size  of  the 
interval.      The  parasitic  solutions  exist. 

Another  important   factor  which   can  cause  numerical  instability 
is  nonlinearity  in  an  ordinary  differential  equation.      We  will  attempt 
to  show  this   in  this  work.      Suppose  that  the   initial  value  problem 
(2.1)  y'    =   f(x,  y)  a  <  x  <  b 

yCa)   =  yQ 

where   f  is   continuous   and  satisfies   Lipschitz  conditions  with 

respect  to  y,   is  to  be   solved  by  a  predictor-corrector  method  of 

the  type  PECE,   i.e.    calculate   a  predicted  value  of  y,   evaluate   f(x,  y) 

at  the  predicted  value,    correct  y,    and  evaluate   f(x,   y)    at  the 

corrected  value.      The   additional  starting  values   of  y  which  are  needed 

to  start   the   solutions   are   assumed  to  have  been   calculated  by  some 

other  method. 

The  predicted  value  yP        of  y   at   x  =   x     ,     is    calculated 
r  "'n+k         J  n+k 

from  the   formula 

(2-2)  Ck  +  «k-l  *n+k-l  +"--+Vn  "  h[bk-l  Wl  +  "-+Vn] 

where  f    stands  for  f(x  ,,,  y   .).y=0,  1,...  k-1.   The  function 
n+y  n+j '  Jn+j   J 

f  is  then  evaluated  at  x  ,,  ,  y  ,,  to  give 

n+k'  Jn+k    & 

fP   =  f(x     vP   ) 
n+k   IV  n+k'  yn+k; 

The  corrected  value  y  ,  of  y  is  then  calculated  from  the  formula. 

Jn+k         J 

(2'3)  ^n+k  +  Ck-1  yn+k-l  +---+C&  =  h[dk  fn+k  +  \-l  W-l  ♦— ♦Vnl- 

Finally,   the   function   f  is    again  evaluated  at   x        ,   y  to  give 

fn+k  =   f(xn+k'   yn+k} 
where  predictor-corrector  pair  is   of  suitable   order. 


Henrici    defines  the  error  in  each  predicted  value  by 

e~  =  y5  -  y(x  )  and  error  in  each  corrected  value  by  e  =  y  -  y(x ) 
m    m      m  m    m      m 

where  m  =  0,  1,  2,...  and  he  shows  that  e~  ,  satisfies 

UM    ** '-  -jio  [Uj  "  hbJ  w  ^j1  "  Vi hP+1  ***ll*J 


*  0(h*+2) 


and  that   en+k  satisfies 


k-1 
(2'5)  en+k  =   ~l      [Cj   "  hdj    Wj]   en+j]   +  hdk  ^n+k  en+k 


"  Cq+1  h4+1  y<l+1(^)   +  °(hq+2) 


In    (2.k)    and   (2.5)   the   a's,   b's,    C's,   and  dfs   are  those   of   (2.2)    and 
(2.3)   and  s's   are   defined  by 

(2.6)  f     _  f(x    ,  y(x   ))   =     5   e     if  e     ^  0 

m  m'  J      m  mm  m 

0     if  e     =0 

m 

The   quantity  C +k  in   (2.5)   is   defined  similarly  by 

(2.7)  f£  -  f(x  ,  y   )   =     5~e~   if  e~  *  0 

m  mm  mm  m 

0     if  e~   =  0 
m 

Also,    in    (2.4)    and   (2.5)   the   constants   C~  . ..    and  C    . .    are  known   for  each 

p+1      q+1 

predictor  and  corrector  method  used.   £~+v  is  eliminated  from  (2.5)  by 
substituting  from  (2.4).  We  get 


-  C     ,   h1+1  y1+1(x  )  +  Odi"1111  (1+2'  p+3)) 
q+1  n 

If  it  is   assumed  that   f(x,  y)  has   continuous   derivatives   of 

suitably  high  order  in  the   interval  of  integration,   then  £'s   of  (2.6) 

are   such  that   r     =   f   (x   ,   t    )  where  t     is  between  y     and  y(x   ).      If 
sm         y     m'     m  m  m  J      m 

t     is   denoted  by  y(x   )   +  Y   ,   then    |y    I    <    le    I    and 
m  ^    ^  x   m  m  m    —       m 

C     =   f   (x    ,  y(x   )   +  v    ) 
^m         y     m'  J      m  'm 

=   t  (x  .  y(x   )   +  Ym  f      (x       y(x   ))   +  O(y^) 
y     m  m  m    yy     mm  m 

Similarly  we   can  show  that 

(2.9)  r   =   f   (x.  y)   +  y~) 

m         y     m     "  m  m 

=  fix,  y(x  J  +  y~  f   r(x  .  y(x  ))  +  O(y^) 
y     m  m  m     yy     mm  m 

where    I v~ I    <    I  e~  I 
1  'm'    —  '    m1 

To  study  the  effect  of  nonlinearity  in  the   differential  equation  on  the 

numerical  stability  of  the  solution  we  proceed  as   follows.      We   define 

(2.10)  g(x)   =   f   (x,   y(x)) 

gm  =   g(xm) 
and  from  (2.9)   and  (2.10) 

(2.11)  r     =   g     +  0(e    ) 

m         m  m 

Since    I  y        <    t  e    I  ,   and 
1    m'    —  ■    m' 


(2.12)  ;m  -  e.  ♦  o(c;) 

Since  |v~l  <    lr~l.   Since  it  is  assumed  that  f(x,  y)  has  continuous 

derivatives  of  suitably  high  order,  it  follows  that 

(2.13)  g^  =  g(xn  +  jh)  =  g^  +  jhg^  +  0(h2)    j  =  0,  1,  2,...,  k-1 

Now,  since  (yn  +  j )  =  0(en  +  j )  =  0(hq)  also  Y~+k  =  °(£n+k)  =  0(hP). 
Now  from  (2.1l)  and  (2.12)  and  (2.13)  we  can  assert  that 
'n+j 


(2.1U)    I  =  gn  +  jh  g£  +  0(h2)    j  =  0,  1,  2,...,  k-1 


***  ^n+j  =  *n  +  *h   *n  +  0(h2) 

In  (2.8)  the  C's  are  replaced  by  the  expressions  given  in  (2.14)  and 
the  terms  rearranged  to  obtain 
k-1 


(2-15)         £n+k  +    4   [lCj  +  h«n  (aA  "  V  "  h \  Vj 


;0  lluj  tngn  iaA-aj'-f2 

J 

h2g;  (Jdj  -  koA'  - h V;  Vj  (j  +  k) 

It    .2 


-  h\  Vi  K+J  +  Vk>  "  ^  Vj  «kVj  +  JW 

-  h6  Vj  Vj  W1  £n+J 

-  c;+i  \  \ +  ^ +  »2  w hP+2  *(p+1)(5) 

-  Cp+1  h*+1  y(1+1)(U   +  0(h^2) 


In  conventional  methods  to  find  the  region  of  stability 
g!    is  taken  to  he   zero.  While  solving  a  differential  equation 

absolute  values   of  g'    and  g     are  not  small  and  thus  may  have   a 

marked  affect  on  the  stability  of  solution.      In  the  literature  one 
can  find  root  locus  plots   of  the   characteristic  equations   of  the 
difference  equation   (2.8)   giving  the  magnitudes   of  roots  plotted 
against  E  *  hg     is   constant.      Krough         gives  these  plots   for  Adams- 

Bashforth  and  Adams -Moulin  type  methods. 

In  this  work  we  will  take   g     to  be   constant.      In  the 

analysis  that  is   given  here,  third  and  fourth  order  terms   in  the 

righthand  side  of   (2.15)    are   disregarded. 

The   characteristic  equation  of  the   simplified  difference 

equation  is 

k— 1 
(2.16)  Yk  +      Z      [[C     +  hgn    (ad     -  d    )    -  hV    (db    ) 


2«I       t  AA         _     V  *      M     ~dl     . 


-h^  ljdj  -kajV]  Y  1  =  0 

_  2 

To  study  the   region  of  stability  we   define  hg     =  h  and  g'/g^   =  ot  in 

(2.16)  and  assume  that   g/    i-  0.      Equation    (2.l6)   thus  becomes 

k       k_1  2 

(2.17)  Y     +     £      [[C,   +  h  (a.d,    -  d.)  -  h"   (d,bj 

i-Q  J  J     -K-  J  K     J 

-  h2   a(jd     -  k    i)]   Y"3]   =  0 

Both  a  and  h  are  taken  to  be  complex  valued.   Of  interest  are  those 

values  of  h  and  a  in  (2.17)  which  will  make  the  roots  y  =   Y->  ip  =  1>  2,...,  k 

of  this  polynomial  satisfy  the  condition  of  stability.   If  |y.  |  £l, 

i  =  1,  2,...,  k  and  roots  of  unit  magnitude  are  simple,  the  method  is 
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said  to  be  absolutely  stable.   The  region  of  absolute  stability  is 

defined  to  be  the  set  of  points  in  h  plane  (for  particular  value  of 

a)  for  which  the  method  is  absolutely  stable. 

One  of  the  roots  '  y  '  of  the  characteristic  equation,  the 

P 

principal  root,   is   approximately  e    .     Any  other  root  we  label  Ye  "to 

indicate  that   it  is   an  extraneous   root.      The  extraneous   roots   result 
from  approximating  a  first  order  differential  equation  with  a 
difference  equation  of  higher  order.      If  each  of  y     satisfies 

|y_|    <  e    ,   and  if  the   y     with  the  magnitude  of  e     are  simple,  then 

the  method  is   said  to  be   relatively  stable.      The   region  of  relative 
stability  is   defined  to  be  the  set  of  points   in  the  h  plane   for 

which  the  method  is   relatively  stable.      Inside  this   region  introduced 

[7] 
errors   grow   (or  decrease)    at  the   same  rate  as  the  solution.      Ralston 

has   defined  a  method  to  be   relatively  stable  provided  all  extraneous 

roots   satisfy    \y        <_  \  y       and  those  y     with  magnitude  of  y     are  simple. 

In  this   region  useful  solutions   can  be  obtained  and  are   computationally 
more   convenient. 

To  keep   computations   at   a  reasonable  level  only  eight  values 
of  a  are   considered  for  each  pair  of  predictor-corrector  methods. 
(2.18)         a  =  0-5eij7r/U         j  =  1,  2,...,  8 

Lehmer  technique  is   used  to  test  whether  all  roots   of  the 

characteristic  equation  lie   inside  the   region    |z|    <,  1.      The  method  is 
as   follows . 
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In  characteristic  equations  substitute  y  =  l/z.   Let 
the  resulting  polynomial  be  denoted  by 

(2.19)  g(z)  =  aQ  +  a1  z  +  ...+8^ 

where  the  coefficients  a.  are  complex. 

Let 

(2.20)  g~(z)  =  a  +  a  +  z  +...+a„  z 

B        n    n  0 

The  linear  combination 

T(g)  =  T(g(z))  =  aQ  g(z)  -  ang  (z) 

has  no  terms  in  z  so  that  T(g)  is  of  lower  degree  than  g.   The 
constant  term  of  T(g)  is  of  particular  interest  because  it  is  real. 
In  fact 


T(g(0))  =  aQa0  -  anan  =  I &r I 


|2      |2 

01 


If  this  constant  term  is  different  from  zero  T(g)  is  a  fit  polynomial 
to  apply  this  transformation.   Continuing  as  far  as  possible  we  obtain 
the  finite  sequence 

T(g),  T2(g),...  -Ag) 
where 

^g(O))  *  0 
If  for  same  i  >  0,  T  (g(0))  <   0,  then  g  has  at  least  one  root  inside 
the  unit  circle.   This  implies  that  equation  (2.17)  has  at  least  one 
root  outside  the  unit  circle.   In  this  way  we  are  able  to  determine 
those  points  for  which  there  is  no  root  of  g(z)  inside  the  unit  circle. 
This  provides  us  the  region  of  absolute  stability.   This  algorithm  also 
provides  a  convenient  method  to  determine  the  region  of  relative 
stability.   If  the  principal  root  =  y     then  we  make  the  transformation 
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z  =   z/y     in  equation   (2.19)   and  follow  the   above  procedure, 


Methods   Investigated 

Adams -Bashforth  predictor  and  Adams-Moulton  corrector  of 
orders  three  through  six  are   considered. 

TABLE  1.      Coefficients   of  Adams -Bashforth  Predictor 

Ck  +  ak-lyiwk-l+---+Vn 

=  h[\-lfn+k-l+---+Vn'  +  VlhP+lyP+1 

a,         =  -1,   all  other  a's   are  zero. 


Fig. 
5-12 

Fig. 
13-20 

Fig. 
21-28 

Fig. 
29-36 
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3 

h 

5 

6 

bo 

5/12 

-9/2U 

251/720 

-U25/1UU0 

bl 

-16/12 

3T/2U 

-127^/720 

2627/lUUo 

b2 

23/12 

-59/24 

2616/720 

-6798/lUUo 

b3 

55/24 

-277^/720 

9U82/1UU0 

\ 

1901/720 

-7673/lUUo 

b5 

U227/1UU0 
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TABLE  2.      Coefficients   of  Adams -Moult on  Corrector 

y         +  C         y  +...+C  y     =  h[d,    f5  ,    +  d,    ,    f     ,    ,    +...+dnf   ] 

°n+k         k-1  yn+k-l  0°n  k     n+k        Tc-1     n+k-J.  On 


♦   C  +_    h*+1  y4+1(xn)    +   0(h*+2) 


.  ,    =  -1,   all  other  C's   are   zero 


Fig. 
5-12 

Fig. 
13-20 

Fig. 
21-28 

Fig. 
29-36 

q. 

3 

1+ 

5 

6 

do 

0 

0 

0 

0 

dl 

-1/12 

1/24 

-19/720 

27/1440 

d2 

2/3 

-5/24 

106/720 

-173/lUUo 

d3 

5/12 

19/24 

-264/720 

482/1440 

\ 

3/8 

646/720 

-798/1440 

d5 

251/720 

1427/1440 

d6 

475/1440 

XX 

NOTE:      x     x     gives  boundary  for  the   region  of  absolute  stability 


•      ■      gives  boundary  for  the   region  of  absolute   and 
relative   stability 


oo 
oo     gives  boundary  for  the  region  of  relative  stability 
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The   curves  of  Figure  1-k  result   from  solving  the   characteristic 
equation   for  a  =  0.      With  the  symmetry  in  mind,   only  the  parts   of  the 
curves  in  the  upper  half  plane   are   given.      Portions  to  the  right  of 
solid  curves   are  the  regions   of  absolute  stability  and  portions  to  the 
right  of  dashed  curves   are  the  regions   of  relative   stability.      For 
a  $  0  stability  was   calculated  at   an  interval  of    «1  along  x-axis   and 
y-axis. 

The   figures   5-36  result   from  solving  the   characteristic 
equation   for  a  =  o»5   ^    '  j   =  1,   2,...    8  for  different 

predictor-corrector  pairs.      Since   in  this    case   there   is   no  symmetry, 
complete   figures   are   given.      When  we   compare  the   region  of  stability 
for  a  =   o   and  for  a  f   0  it   is   quite   obvious  that   even   a  small  value 
of  a  has   a  marked  effect   on  the  region  of  stability • 

For  any  predictor-corrector  pair  it  was    conjectured  that 
regions   of  absolute    (relative)   stability  for  a  #  0  will   always   lie 
inside  the   region   for  a  =  0.      But  this   is   not  true.      Compare  Figure  22 
with  Figure   3.      In  Figure   3  the  point  0  +    . 8i   is  not   stable  where   as   in 
Figure  22  point  0  +    . 8i   is   stable.      But  in  general,   the   region  of 
stability  shrinks  when  a  4-  0. 

Another  conjecture  which  is   disproved  is  the   following.      We 
have   defined  a  function 

6(a,  h)    s0  if  stable 

1  if  unstable 
6(0,   5(1  -   jah))   =0  j   =   0,   1,...    k 

=>      9(a,   h)   =  0 
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Absolute  Stability.   In  Figure  lU  o  a  ..'5i.   The  point 
E  =  -.1  +  .3i  is  unstable  or  8(a,  E)  =  1 

but  6(0,  5(1  -  Jah))  =0    j  =  0,  1,...  h 

from  Figure  2.  k  =  h   for  both  of  these  figures. 

Relative  Stability.   Consider  again  Figure  lU,  a  =  .5i. 
The  point  E  =  +.1  -  . i  is  unstable 

but  6(0,  E(l  -  jotE))  =  0    J  =  0,  1,...  k 

from  Figure  2. 

Another  interesting  fact  to  note  from  the  study  of  these 
figures  is  as  the  order  of  method  is  increased  the  region  of  stability 
narrows  down  considerably. 


3k 
3.      SUMMARY 

For  any  specified  predictor-corrector  pair  the   a's,  b's,   C's, 
and  d's   in   (2.17)   are  known.      We  have  seen  that   stability  depends   upon 
the   roots   of  characteristic  equations.     We  have  replaced 
«      .   =   g(x     +j)   =  £L+  j^gl.     We  neglect  the  higher  order  terms.      It 

is  not  clear  how  the  neglected  terms  will  effect  the  situation.     When 
we   compare  the  region  of  stability  for  a  ^  0  with  a  =  0  it   is   evident 
that  the   region  of  stability  shrinks   in  most   cases.      Suppose  that   in 
the   course  of  the  solution,   estimates   of  g     and  g*    for  each  n  are 

known  and  it   is   found  that  neither  is   small  in  absolute  value.      Then 
it  is  very  likely  that  the  variation  in  g     and  g*    can  have   a  marked 

effect  on  the  numerical  stability  of  the  solution.      In  higher  order 
methods  the   region  of  stability  is   quite  smaller  than  the  region  of 
stability  for  lower  order  methods.      Through  the  use  of  Figures   U-36, 
it   is  possible  to  predict  stability  characteristics   of  predictor- 
corrector  methods   as   a  function  of  h  and  a. 
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